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Pebble-bed nuclear reactor technology, which is currently being revived around the world, raises 
fundamental questions about dense granular flow in silos. A typical reactor core is composed of 
graphite fuel pebbles, which drain very slowly in a continuous refueling process. Pebble flow is 
poorly understood and not easily accessible to experiments, and yet it has a major impact on 
reactor physics. To address this problem, we perform full-scale, discrete-element simulations in 
realistic geometries, with up to 440,000 frictional, viscoelastic 6cm-diameter spheres draining in a 
cylindrical vessel of diameter 3.5m and height 10m with bottom funnels angled at 30° or 60°. We also 
simulate a bidisperse core with a dynamic central column of smaller graphite moderator pebbles and 
show that little mixing occurs down to a 1:2 diameter ratio. We analyze the mean velocity, diffusion 
and mixing, local ordering and porosity (from Voronoi volumes), the residence-time distribution, 
and the effects of wall friction and discuss implications for reactor design and the basic physics of 
granular flow. 



I. INTRODUCTION 



A. Background 



A worldwide effort is underway to develop more eco- 
nomical, efficient, proliferation resistant, and safer nu- 
clear power f3|. A promising Generation IV reactor de- 
sign is the uranium-based, graphite moderated, helium- 
cooled very high temperature reactor 0, which of- 
fers meltdown-proof passive safety, convenient long-term 
waste storage, modular construction, and a means of 
nuclear-assisted hydrogen production and desalination. 
In one embodiment, uranium dioxide is contained in mi- 
crospheres dispersed in spherical graphite pebbles, the 
size of billiard balls, which are veryslowly cycled through 
the core in a dense granular flow 0| . Control rods are 
inserted in graphite bricks of the core vessel, so there are 
no obstacles to pebble flow. 

The pebble-bed reactor (PBR) concept, which orig- 
inated in Germany in the 1950s, is being revisited by 
several countries, notably China "F] (HTR-10 ''^) and 
South Africa 3] (PBMR which plan large-scale de- 
ployment. In the United States, the Modular Pebble 
Bed Reactor (MPBR) jllsl is a candidate for the Next 
Generation Nuclear Plant of the Department of Energy. 
A notable feature of MPBR (also present in the origi- 
nal South African design) is the introduction of graphite 
moderator pebbles, identical to the fuel pebbles but with- 
out the uranium microspheres. The moderator pebbles 
form a dynamic central column, which serves to flatten 
the neutron flux across the annular fuel region without 
placing any fixed structures inside the core vessel. The 
annular fuel region increases the power output and effi- 
ciency, while preserving passive safety. In the bidisperse 
MPBR, the moderator pebbles are smaller to reduce the 
permeability of the central column and thus focus helium 
gas on the outer fuel annulus. The continuous refueling 
process is a major advantage of pebble-bed reactors over 



other core designs, which typically require shutting down 
for a costly dismantling and reconstruction. The random 
cycling of pebbles through a flowing core also greatly im- 
proves the uniformity of fuel burnup. 

In spite of these advantages, however, the dynamic 
core of a PBR is also a cause for concern among de- 
signers and regulators, since the basic physics of dense 
granular flow is not fully understood. Indeed, no reliable 
continuum model is available to predict the mean veloc- 
ity in silos of different shapes 0, although the empirical 
Kinematic Model 

ES in 111 provides a reasonable fit 
near the orifice in a wide silo 0, lisl IT^ . A mi- 
croscopic model for random-packing dynamics has also 
been proposed 0| and fitted to reproduce drainage in a 
wide silo 0, but a complete statistical theory of dense 
granular flow is still lacking. The classical kinetic theory 
of gases has been successfully applied to dilute granular 
flows I23, 01 , in spite of problems with inelastic colli- 
sions '22*1 , but it clearly breaks down in dense flows with 
long-lasting, frictional contacts |16, 23], as in pebble-bed 
reactors. Plasticity theories from soil mechanics might 
seem more appropriate jl2l |. but they cannot describe 
flows in silos of arbitrary shape and often lead to violent 
instabilities p4L l25l | , although a stochastic flow rule |2^ 
may resolve these difficulties and eventually lead to a 
general theory. 

For now, experiments provide important, although lim- 
ited, information about dense granular flows. Many ex- 
periments have been done on drainage flows in quasi-2d 
silos where particles are tracked accurately at a transpar- 
ent wall I3, EJ, Ea, Ea I23 • Some three-dimensional par- 
ticle tracking in granular materials and colloids has also 
been done with magnetic resonance imaging confo- 
cal microscopy j22j, index matching with an interstitial 
fluid , and diffusing- wave spectroscopy ^sT] , although 
these systems are quite different from a pebble-bed re- 
actor core. Experimental studies of more realistic ge- 
ometries for PBR have mostly focused on the porosity 
distribution of static packings of spheres |33. l3^ , which 
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aflFects helium gas flow through the core |3 HE HBl ■ 

As a first attempt to observe pebble dynamics experi- 
mentally in a reactor model, the slow flow of plastic beads 
has recently been studied in 1:10 scale models of MPBR 
in two different ways |^] : The trajectories of colored peb- 
bles were recorded (by hand) along a plexiglass wall in 
a half-core model, and a single radioactive tracer pebble 
in the bulk was tracked in three dimensions in a full-core 
model. Very slow flow was achieved using a screw mech- 
anism at the orifice to approximate the mean exit rate 
of one pebble per minute in MPBR. These experiments 
demonstrate the feasibility of the dynamic central column 
and confirm that pebbles diffuse less than one diameter 
away from streamlines of the mean flow. However, it is 
important to gain a more detailed understanding of peb- 
ble flow in the entire core to reliably predict reactor power 
output, fuel efficiency, power peaking, accident scenarios 
using existing nuclear engineering codes js^, l39l | . 



B. Discrete-Element Simulations 

Simulations are ideally suited to provide complete, 
three-dimensional information in a granular flow. Some 
simulations of the static random packin g of fuel pebbles 
in a PBR core have been reported j40l l4lj. but in the 
last few years, large-scale, parallel computing technol- 
ogy has advanced to the stage where it is now possi- 
ble to carry out simulations of continuous pebble flow 
in a full-sized reactor geometry using the Discrete Ele- 
ment Method (DEM). In such simulations, each particle 
is accurately modeled as a sphere undergoing realistic 
frictional interactions with other particles ^3 In 
this paper, we present DEM simulations which address 
various outstanding issues in reactor design, such as the 
sharpness of the interface between fuel and moderator 
pebbles (in both monodisperse and bidisperse cores), the 
horizontal diffusion of the pebbles, the geometry depen- 
dence of the mean streamlines, the porosity distribution, 
wall effects, and residence-time distributions. 

Our simulations are based on the MPBR geome- 
try 0, , consisting of spherical pebbles with diameter 
d = 6cm in a cylindrical container approximately 10m 
high and 3.5m across. In this design there is a central 
column of moderating reflector pebbles, surrounded by 
an annulus of fuel pebbles. The two pebble types are 
physically identical except that the fuel pebbles contain 
sand-sized uranium fuel particles. Particles are continu- 
ously cycled, so that those exiting the container are rein- 
troduced at the top of the packing. In order to efficiently 
maintain the central column, a cylindrical guide ring of 
radius rjn = 14.5c? extends into the packing to z = 140(i. 
Reflector pebbles are poured inside, while fuel pebbles 
are poured outside, and the guide ring ensures that two 
types do not mix together at the surface. Figure^shows 
the two main geometries that were considered; for much 
of this analysis, we have concentrated on the case when 
the exit funnel is sloped at thirty degrees, but since this 



angle can have a large effect on the pebble flow, we also 
consider the case of the when the funnel is sloped at sixty 
degrees. In both cases the radius of the opening at the 
bottom of the funnel is roxit = 5(i. 

In MPBR, as in most pebble-bed reactors, the drainage 
process takes place extremely slowly. Pebbles are indi- 
vidually removed from the base of the reactor using a 
screw mechanism, at a typical rate of one pebble per 
minute, and the mean residence time of a pebble is 77 
days. Carrying out a DEM simulation at this flow rate 
would make it infeasible to collect enough meaningful 
data. However, previous experimental work by Choi et 
al. [l^] has shown that the regime of slow, dense granu- 
lar flow is governed by a distinctly non-thermal picture, 
where particles undergo long-lasting contacts with their 
neighbors, and the features of the flow are predominately 
governed by geometry and packing constraints. In par- 
ticular, they observed that for a large range of hopper 
drainage experiments, altering the orifice size resulted in 
a change in the overall flow rate, but did not alter the ge- 
ometry of the flow proflle - the flow velocities were scaled 
by a constant factor. Furthermore, geometric properties 
of the flow, such as particle diffusion, were unaffected 
by the overall flow rate. We therefore chose to study a 
faster flow regime in which pebbles drain from the reac- 
tor exit pipe under gravity. Our results can be related 
directly to the reactor design by rescaling the time by an 
appropriate factor. 

As well as the two full-scale simulations described 
above, we also considered a half-size geometry in order 
to investigate how various alterations in the makeup of 
the reactor would affect the flow. In particular, we exam- 
ined a series of bidisperse simulations, in which the di- 
ameter of moderator particles in the central column was 
reduced. As explained in section IVIIII this has the effect 
of reducing the gas permeability of the central column, 
thus focusing the helium coolant flow on the hottest re- 
gion of the reactor core, in and around the fuel annulus. 
The purpose of the simulations is to test the feasibility 
of the bidisperse PBR concept, as a function of the size 
ratio of moderator and fuel pebbles, with regard to the 
granular flow. It is not clear a priori under what condi- 
tions the dynamic column will remain stable with little 
interdiffusion of moderator and graphite pebbles. 

To study this issue, we made a sequence of three runs 
using a half-size reactor geometry. (The smaller core size 
is needed since the number of smaller pebbles increases 
as the inverse cube of the diameter ratio.) The geometry 
is similar to that used above, except that the radius of 
the cylindrical container is decreased to 15(i, with the 
guide ring at r„^ — 7.5d. The radius of the exit pipe 
is decreased to roxit = 4,d. In the experiments, we keep 
the diameter of the fuel pebbles flxed at d, and use d, 
0.8d, and 0.5d for the diameters of the moderator pebbles. 
The same geometry was also used to study the effect 
of wall friction, by making an additional run with the 
particle/wall friction coefficient /i^ = 0. 

The paper is organized as follows. In section^ we dis- 
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FIG. 1: (Color online) Snapshots of vertical cross-sections of 
the simulations for the two geometries considered in this re- 
port. We make use of a cylindrical coordinate system (r, 9, z) 
where 2 = at the orifice. At the base of the container there is 
a small exit pipe of radius rexit = 5ci that extends upwards to 
z = lOd. This connects to a conical funnel region, which has 
slope thirty degrees (left) or sixty degrees (right). The con- 
ical wall connects to a cylindrical wall of radius rout ~ 29d, 
a,t z = 23.86ci and z = 51.57d for the thirty and sixty degree 
reactor geometries respectively. Particles are poured into the 
container up to a height of approximately z = 160d. A cylin- 
drical wall at rin — 14. 5d extends down into the packing to a 
height oi z = 140d to keep the two types of pebbles mixing 
at the surface. 



CUSS the simulation technique that was used and briefly 
describe its implementation. This is followed with some 
basic analysis of the velocity profiles and a comparison to 
the Kinematic Model in section Hill We study diffusion 
around streamlines in section HV1 and the distribution of 
porosity and local ordering in section^ Next, in section 
I VII we examine the residence-time distribution of pebbles 
in the reactor, which is related to fuel burnup, and in 
section I VIII we show that wall friction plays an impor- 
tant role. In section IVIIII we analyze the bidisperse PBR 
concept with half-size reactor simulations for a range of 
pebble-diameter ratios, focusing on the mean flow, dif- 
fusion, and mixing. We conclude in section Hxl by sum- 
marizing implications of our study for reactor design and 
the basic physics of granular flow. 



II. MODELS AND METHODS 

The DEM simulations are based on a modified version 
of the model developed by Cundall and Strack to 
model cohesionless particulates 0|. Monodisperse 
spheres with diameter d interact according to Hertzian, 
history dependent contact forces. For a distance r be- 
tween a particle and its neighbor, when the particles are 
in compression, so that (5 = d — |r| > 0, then the two par- 
ticles experience a force F = F„ -|~ Ft , where the normal 
and tangential components are given by 

F„ = ^d(krM-^) (1) 

Ft = VVrf (-fc* Ast - ^) . (2) 

Here, n — r/|r|. v„ and vt are the normal and tan- 
gential components of the relative surface velocity, and 
kn,t and ^n,t are the elastic and viscoelastic constants, 
respectively. Ast is the elastic tangential displacement 
between spheres, obtained by integrating tangential rela- 
tive velocities during elastic deformation for the lifetime 
of the contact, and is truncated as necessary to satisfy a 
local Coulomb yield criterion |Ft| < /i |F„|. Particle- wall 
interactions are treated identically, though the particle- 
wall friction coefficient /z„ is set independently. 

For the monodispersed system, the spheres have diam- 
eter d = 6cm, mass m = 210g and interparticle friction 
coefficient /i = 0.7, flowing under the influence of gravity 
g = 9.81ms~^. For the bi-dispersed systems, the modera- 
tor particles have diameter 0.8c? or O.Sd. The particle- wall 
friction coefficient /i^ = 0.7 except in one case where we 
model a frictionless wall, jiw = 0.0. For the current sim- 
ulations we set fct = |A;„, and choose fc„ = 2 x lO^mg/d. 
While this is significantly less than would be realistic 
for graphite pebbles, where we expect fc„ > lO^^mg/d, 
such a spring constant would be prohibitively computa- 

— 1/2 

tionally expensive, as the time step scales as 6t 
for collisions to be modeled effectively. Previous simu- 
lations have shown that increasing fc„ does not signifi- 
cantly alter physical results '4^. We use a time step of 
5t = 1.0 X 10~'*r and dam ping coefficients 7„ = 50t~^ 
and 7t = 0.0, where r = \fdfg = 0.078s. All measure- 
ments are expressed in terms of d, m and r. 

The initial configurations are made by extending the 
inner cylinder from 140c? to the bottom of the container, 
adding a wall at the bottom of the container to stop par- 
ticles from draining, and pouring in moderator pebbles 
into the inner cylinder and fuel pebbles between the in- 
ner and outer cylinders until the reactor was loaded. The 
bottom wall is then removed, the inner cylinder is raised 
to 140c?, and particles are allowed to drain out of the 
container. As noted above, particles are recycled with 
moderator particles reinserted within the inner cylinder, 
and fuel particles between the inner and outer cylinders. 
All results presented here are after all the particles have 
cycled through the reactor at least once. The number of 
moderator and fuel particles was adjusted slightly from 
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the initial filling so that the level at the top of the reac- 
tor is approximately equal. For the full scale simulation 
with a thirty degree outlet, the total number of pebbles 
is 440,000 with 105,011 moderator pebbles and 334,989 
fuel pebbles, while for the sixty degree outlet, the to- 
tal number of pebbles is 406,405 with 97,463 moderator 
and 308,942 fuel pebbles. For the former case, a million 
steps took approximately 13 hours on 60 processors on 
Sandia's Intel Xenon cluster. 

For the bidispersed simulations the total number of 
pebbles is 130,044, 160,423, and 337,715 for the diame- 
ter of the moderator particles equal to d, 0.8d and 0.5d 
respectively. As the diameter of the moderator pebbles 
is decreased the number of particles required rapidly in- 
creases, since it scales according to the inverse of the 
diameter cubed. 

A snapshot of all the particle positions is recorded ev- 
ery 5r = 0.39s. For the thirty degree reactor geometry we 
collected 1,087 successive snapshots, totaling 24.9Gb of 
data, while for the sixty degree reactor geometry, we col- 
lected 881 successive snapshots, totaling 18.7Gb of data. 
A variety of analysis codes written in Perl and G-I--I- were 
used to sequentially parse the snapshot files to investigate 
different aspects of the flow. We also created extended 
data sets, with an additional 440 snapshots for the thirty 
degree geometry, and 368 snapshots for the sixty degree 
geometry, for examining long residence times in section 



III. MEAN- VELOCITY PROFILES 
A. Simulation Results 

Since we have a massive amount of precise data about 
the positions of the pebbles, it is possible to reconstruct 
the mean flow in the reactor with great accuracy. How- 
ever care must be taken when calculating velocity profiles 
to ensure the highest accuracy. Initial studies of the data 
showed that crystallization effects near the wall can cre- 
ate features in the velocity profile at a sub-particle level, 
and we therefore chose a method that could resolve this. 

By exploiting the axial symmetry of the system, one 
only need to find the velocity profile as a function of 
r and z only. The container is divided into bins and 
the mean velocity is determined within each. A particle 
which is at x„ at the nth timestep and at x„_|_i at the {n+ 
l)th timestep, makes a velocity contribution of (x„_|_i — 
x„)/Ai in the bin which contains its midpoint, (x„+i + 
x„)/2. 

In the z direction, we divide the container into strips 
Id across. However, in the r direction we take an alter- 
native approach. Since the number of pebbles between a 
radius of r and r -I- Ar is proportional to rAr, dividing 
the container into bins of a fixed width is unsatisfactory, 
since the amount of data in bins with high r would be 
disproportionately large. We therefore introduce a new 
coordinate s = r^. The coordinate s covers the range 
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FIG. 2: (Color online) Computed streamlines of the mean 
flow in the 30° (left) and 60° reactor geometries. Arrows 
are proportional to the velocity vectors in selected horizontal 
slices. 
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and we divide the container into regions 



that are equally spaced in s, of width Id^. The number 
of pebbles in each bin is therefore roughly equal, allowing 
for accurate averaging in the bulk and high resolution at 
the boundary. 

This result yields extremely accurate velocity profiles 
in the cylindrical region of the tank. However, it fails to 
capture crystallization effects in the conical region: since 
the particles are aligned with the slope of the walls are 
averaged over a strip in z of width Id, any effects are 
smeared out across several bins. We therefore scaled the 
radial coordinate to what it would be if the particle was 
in the center of the strip. Specifically, if the radius of 
the container is given by R{z), a particle at (r„,z„) is 
recorded as having radial coordinate rnR{z) / R{zn). In 
the cylindrical region of the tank this has no effect, while 
in the conical region, it effectively creates trapezoid- 
shaped bins from which it is easy to see crystallization 
effects which are aligned with the wall. 

The streamlines of the mean flow are shown in Fig. [21 
in the two geometries. Streamlines are computed by La- 
grangian integration of the DEM velocity field, starting 
from points at a given height, equally spaced in radius. 
In each geometry, there is a transition from a nonuni- 
form converging fiow in the lower funnel region to a 
nearly uniform plug flow in the upper cylindrical region, 
consistent with the standard engineering picture of silo 
drainage In the wider funnel, there is a region of 

much slower flow near the sharp corner at the upper edge 
of the funnel. Our results for both geometries are quite 
consistent with particle-tracking data for quasi-2d silos of 
similar shapes |3] and half-cylinder models of the MPBR 
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FIG. 3: (Color online) Velocity profiles for the thirty degree 
reactor geometry for several low cross-sections (a) and several 
high cross-sections (b). 



core |37| . which provides an important vaUdation of our 
simulations. 

We now look more closely at horizontal slices of the 
velocity field. Figure |3Ja) shows several velocity profiles 
for the thirty degree case in the narrowing section of the 
container. As expected, we see a widening of the velocity 
profile as z increases. We can also see lattice effects, 
spaced at \Jid apart, due to to particles crystallizing on 
the conical wall section. 

Figure |2Ib) shows similar plots for several heights in 
the upper region of the container. At these heights, the 
velocity profile is roughly uniform across the container. 
However a boundary layer of slower velocities, several 
particle diameters wide, still persists. The average veloc- 
ities of particles touching the boundary is between one 
half and two thirds that of particles in the bulk; it is ex- 
pected that this behavior is very dependent on particle- 
wall friction; this issue is studied in more detail in section 

rviTi 



High in the container, results for the sixty degree ge- 
ometry are very similar to the thirty degree case (and 
thus are not shown). However, as would be expected, 
a significantly different crossover from parabolic fiow to 
plug-like flow in the lower part of the tank is observed, 
as shown in figure [S] 
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FIG. 4; (Color online) Streamlines of the mean flow in the 30° 
(left) and 60° reactor geometries for the numerical solution of 
the Kinematic Model. Arrows are proportional to the velocity 
vectors in selected horizontal slices. 



B. Comparison with the Kinematic Model 

Perhaps the only continuum theory available for the 
mean flow profile in a slowly draining silo is the Kine- 
matic Model [TSinillHill, which postulates that hori- 
zontal velocity vector u is proportional to the horizontal 
gradient V_l of the downward velocity component v (i.e. 
the local shear rate), 
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where h is the "diffusion length" , a material parameter 
typically in the range of one to three particle diameters. 
The idea behind Eq. ^ is that particles drift from re- 
gions of low to high downward velocity, where there are 
more local rearrangements (and more free volume) to ac- 
commodate their collective motion. The approximation 
of incompressibility, V • (u, —v) = 0, applied to Eq. (PJ 
yields a diffusion equation for the downward velocity. 



(4) 



where the vertical coordinate z acts like "time" . Bound- 
ary conditions on Eq. I^J require no normal velocity com- 
ponent at the container walls, except at the orifice, where 
V is specified (effectively an "initial condition"). As de- 
scribed in Appendix^ this boundary-value problem can 
be accurately solved using a standard Crank-Nicholson 
scheme for the diffusion equation. 

The kinematic parameter 6 can be understood as a dif- 
fusion length for free volume, which is introduced at the 
orifice and diffuses upward, causing downward diffusion 
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FIG. 5: (Color online) Velocity profiles for the 60° reactor 
geometry (solid lines), with a comparison to the Kinematic 
Model for b — 3d (dashed lines). 



to capture the rapid transition from converging flow to 
plug flow seen in the DEM data. This is shown clearly 
by comparing the streamlines for the Kinematic Model 
in figure0]with those for DEM. Streamlines for the Kine- 
matic Model are roughly parabolic, and no single value of 
b can capture the rapid change from downward streani- 
hnes to converging streamlines seen in DEM. 

The difflculty in precisely determining b is also a com- 
mon theme in experiments, although recent data suggests 
that a nonhnear diffusion length may improve the fit ^ . 
Perhaps a more fundamental problem with the Kinematic 
Model is that it cannot easily describe the rapid crossover 
from parabolic converging flow to uniform plug flow seen 
in both geometries our DEM simulations; we will return 
to this issue in section Ivl 



IV. DIFFUSION AND MIXING 



of particles. It was originally proposed that free volume is 
carried by voids fiol llsf , which displace single particles 
as they move, but a more realistic mechanism involves 
cooperative particle motion due to diffusing "spots" of 
delocalized free volume [13. The Spot Model can pro- 
duce accurate flowing packings in wide silos |0| , and the 
Kinematic Model can be derived as the continuum limit 
of the simplest case where spots drift upward at constant 
velocity (due to gravity) while undergoing independent 
random walks, although more general continuum equa- 
tions are also possible for different spot dynamics. A 
first-principles mechanical theory of spot dynamics is still 
lacking (although it may be based on a stochastic refor- 
mulation of Mohr- Coulomb plasticity [2^), so here we 
simply try a range of b values and compare to the DEM 
flow proflles. 

Consistent with a recent experimental study of quasi- 
2d silos 1^, we find reasonable agreement between the 
Kinematic Model predictions and the DEM fiow pro- 
files, but the effect of the container geometry is not 
fully captured. In the converging flow of the funnel 
region, the streamlines are roughly parabolic, as pre- 
dicted by the Kinematic Model and found in many exper- 
iments OE^QjEiEI- For that region, it is possible to 
choose a single value (6 = 3d) to achieve an acceptable fit 
to the DEM fiow profiles for both the 30° and 60° funnel 
geometries, as shown in figure |31 

In spite of the reasonable overall fit, the Kinematic 
Model has some problems describing the DEM results. It 
fails to describe the several particle thick boundary layer 
of slower velocities seen in the DEM data. In the original 
model, b depends only on the properties of the granular 
material, but we find that it seems to depend on the ge- 
ometry; the best fit to the 30° DEM data is 6 « 2.5d, 
while the best fit for the 60° DEM data is 6 « 3.Qd. Such 
discrepancies may partly be due to the boundary layers, 
since in the lower section of the container the conical 
walls may have an appreciable effect on the majority of 
the flow. We also flnd that the Kinematic Model fails 



Nuclear engineering codes for PBR core neutronics 
typically assume that pebbles flow in a smooth laminar 
manner along streamlines, with very little lateral diffu- 
sion [s^ [s^. Were such significant diffusion to occur 
across streamlines, it could alter the core composition 
in unexpected ways. In the MPBR design with a dy- 
namic central column Q , diffusion leads to the unwanted 
mixing of graphite pebbles from the central reflector col- 
umn with fuel pebbles from the outer annulus, so it must 
be quantifled. Simulations and experiments are crucial, 
since diffusion in slow, dense granular flows is not fully 
understood 17]. 

Particle-tracking experiments on quasi-2d silos [l6j | and 
half-cylinder MPBR models [33 have demonstrated very 
little pebble diffusion in slow, dense flows, but the obser- 
vations were made near transparent walls, which could 
affect the fiow, e.g. due to ordering (see below). Three- 
dimensional tracking of a radioactive tracer in a cylindri- 
cal MPBR model has also shown very little diffusion, at 
the scale of a single pebble diameter for the duration of 
the flow [33. Here, we take advantage of the complete 
information on pebble positions in our DEM simulations 
to study core diffusion and mixing with great accuracy. 

We collected extensive statistics on how much pebbles 
deviate from the mean-flow streamlines during drainage. 
Consistent with theoretical concepts [l^, experiments 
have demonstrated that the dynamics are strongly gov- 
erned by the packing geometry, so that diffusion can most 
accurately be described by looking at the mean-squared 
horizontal displacement away from the streamline, as a 
function of the distance dropped by the pebble (not time, 
as in molecular diffusion), regardless of the flow rate. Mo- 
tivated by the importance of quantifying mixing at the 
fuel/moderator interface in the dynamic central column 
of MPBR, we focus on tracking pebbles passing through 
z — llOd with |r — 15d| < O.lQd. The variance of the r 
coordinate of the particles as they fall to different heights 
in z can be calculated. From this, we can determine the 
amount of radial diffusion, defined as the increase in the 
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V. PACKING STATISTICS 
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FIG. 6: (Color online) Radial diffusion of particles about 
streamlines of the mean flow as a function of height, z, in 
both reactor geometries for pebbles starting ai z — llOd in 
an annulus of radius r — 15d, at the edge of the dynamic 
central column in MPBR. For the 30° geometry, we also show 
data for pebbles near the wall at r = 25d. 



variance of r of the tracked particles from the variance 
at the initial height. 

The diffusion data for both reactor geometries is shown 
in figure IHl We see that for large values of z in the cylin- 
drical part of the container, the pebbles undergo essen- 
tially no diffusion; this is to be expected, since we have 
seen that in this area the packing is essentially plug-like, 
and particles are locked in position with their neighbors. 
However for lower values of z the amount of radial spread- 
ing begins to increase, as the particles experience some 
rearrangement in the region corresponding to converging 
flow. Note however that the scale of this mixing is very 
small, and is much less than a pebble diameter. For very 
small values of z, there is a decrease in the variance of 
the radial coordinate, since the pebbles must converge on 
the orifice as they exit the container. 

We apphed a similar analysis for different initial val- 
ues of r, and found very similar results over the range 
< r < 25d. However, for particles close to the con- 
tainer boundary, very different behavior is observed, as 
shown by the third line in figure |^ for particles with 
\r — 25d\ < O.lOd. In this region, the particles undergo 
rearrangement, and this causes a (piecewise) linear in- 
crease in the mean-squared displacement with distance 
dropped, which corresponds to a constant local diffusion 
length. There is also evidence of a sharp transition in the 
boundary-layer diffusion length, which increases signifi- 
cantly as pebbles pass the corner into the converging-flow 
region of the funnel. 



A. Pebble Volume Fraction 

Pebble-bed experiments [s^, Is^l and simulations 
of static sphere packings in cylinders have revealed 
that there are local variations in porosity near walls, at 
the scale of several pebble diameters, but there has been 
no such study of flowing packings, averaging over dy- 
namic configurations. Similar findings would have impor- 
tant implications for helium flow in the core, since the lo- 
cal gas permeability is related to the porosity [33.l35ll3^ . 

First, we study the distribution of local volume frac- 
tion (% of volume occupied by pebbles) throughout the 
container, averaged in time. (The porosity is one minus 
the volume fraction.) Random close packing of spheres 
corresponds to a volume fraction in the range 55% - 63%, 
while flows occur in a somewhat more restricted range. 
The lower bound is approximately set by random loose 
packing, where rigidity percolation sets |46|, while the 
upper bound is near the jam ming point [47| or the max- 
imally random jammed state where flow cannot oc- 
cur. 

The best way to determine the volume fraction on 
a local scale is to use a Voronoi tessellation, which we 
compute with a novel efficient algorithm for flows, to 
be described in detail elsewhere. The Voronoi tessella- 
tion uniquely assigns a polygonal volume to each pebble, 
formed by intersecting the planes bisecting the lines be- 
tween different pebble centers. The local packing fraction 
in a small region can then be found by taking the ratio 
of the particle volume in that region to the ratio of the 
Voronoi volume. Such a method can be used to define 
local density even down to the scale of a single particle, 
but for this work we compute local density by averaging 
on a scale of several particle diameters. 

Figure [7| shows density snapshots for cross sections 
through the thirty degree and sixty degree reactor ge- 
ometries, based on computing the local density at a par- 
ticular point by averaging over the Voronoi densities of 
particles within a radius of 2.2d. Figure |S1 shows density 
plots over the entire flow of the data, but using a smaller 
averaging radius of 0.8d. Many interesting features are 
visible, which corroborate our other results. High in the 
center of the container, we see that the local packing frac- 
tion is mostly close to 63%, suggesting that the plug-like 
region is in a nearly jammed and rigid state. This is con- 
sistent with our earlier data showing nearly uniform plug 
flow with no significant diffusion or mixing. 

We also observe two annular lines of lower density 
propagating down from the guide ring, which form due 
to wall effects on the guide ring itself (see below) and are 
advected downward. The fact that these subtle artifacts 
of the guide-ring constraints are felt far down in the flow 
further demonstrates that very little diffusion or shear- 
ing occurs in the upper region. There are also similar 
lower-density regions along the walls, related to partial 
crystallization described in more detail below. 




FIG. 7: (Color online) Plots of local volume fraction (1 — 
porosity) in a vertical cross section for the thirty degree re- 
actor geometry (left) and the sixty degree reactor geometry 
(right), calculated using a Voronoi cell method. The color 
scheme used is red 50%, yellow 57%, dark blue 60%, cyan 
63%. High in the bulk of the container, the packing fraction 
is approximately 63%, apart from in a small region of lower 
density at rin = 14. 5d, corresponding to packing defects intro- 
duced by the guide ring. In both geometries a sharp reduction 
in density is observed in a region above the orifice, where par- 
ticles in the parabolic fiow region are forced to undergo local 
rearrangements. 



It is also clear in both geometries, especially the 30° 
model, that there is a fairly sharp transition between the 
upper region of nearly rigid plug flow and a less dense 
lower region of shear flow in the funnel. Similar features 
are in the velocity profiles described above, but the tran- 
sition is much more sharp, at the scale of at most a few 
particles, for the local packing fraction. These sudden 
variations in material properties and velocities are rem- 
iniscent of shock-like discontinuities in Mohr-Coulomb 
plasticity theories of granular materials 0, 0] . It seems 
no such existing theory can be applied to the reactor 
flows, but our results suggest that plasticity concepts 
may be useful in developing a continuum theory of dense 
granular flow p^ . 
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FIG. 8: (Color online) Time- averaged plots of the local vol- 
ume fraction, using the same color scheme as figure |7| 



B. Local Ordering and Porosity 

As noted above, previous simulation studies of local or- 
dering near walls have focused on static packings in sim- 
plifled cylindrical geometries (without the funnel, outlet 
pipe, or guide ring) [iol l4ll|. while we compute average 
statistics for slowly flowing packings in realistic full-scale 
reactor models. To take a closer look at ordering near 
walls, we study the number density profile in horizon- 
tal slices at different heights. The container is divided 
into bins in the same way as discussed previously and 
the number density in a bin is obtained by counting the 
number of times a particle center lies within that bin. 

Figure El^a) shows a sequence of number density pro- 
files for several low values of z in the thirty degree re- 
actor geometry. At all four heights, lattice effects are 
clearly visible and quite similar to those observed in ex- 
periments [S^, nil and other simulations [40l I4H . For the 
lowest three heights, these peaks are roughly -s/Sd apart, 
corresponding to particles crystallized against the coni- 
cal wall, while for the highest value of z, these effects 
are roughly Id apart, due to particles being crystallized 
against the cylindrical wall. The above graph also shows 
that in the middle of the container, no lattice effects are 
present. 

However, this situation changes dramatically higher up 
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FIG. 9; (Color online) Number density plots in the thirty 
degree reactor geometry for several low cross sections (a), and 
several high cross sections (b). 
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of the cylindrical part of the reactor. At much lower 
heights, around z — 40d, this radial ordering is broken, 
as the particles are forced to reorganize once they enter 
the parabolic region of flow. 

To make a direct connection with the modeling of gas 
flow, we show horizontal slices of the porosity at different 
heights in figure 1101 The porosity is measured here by 
intersecting the spheres with annular cylindrical bins to 
compute the fraction of each bin volume not occupied by 
pebbles. The features noted above appear in the porosity 
and alter the local permeability, which enters continuum 
descriptions of helium gas flow in the core |34i fSSl. l36| . 



VI. RESIDENCE-TIME DISTRIBUTION 
A. Predictions of the Kinematic Model 

The statistical distribution of fuel burnup is closely re- 
lated to the distribution of pebble residence times in the 
reactor core, differing only due to nonuniform sampling 
of the neutron flux proflle. Since the upper pebble flow 
is essentially a uniform plug flow, the distribution of res- 
idence times is the same (up to a constant time shift) 
as the distribution of waiting times for pebbles starting 
at a given height in the core to exit through the orifice, 
and we concentrate on these distributions in this section. 
However, we conclude by examining the residence times 
for particles to pass through the entire container, to in- 
vestigate the effects of the guide ring and the outer walls. 

We have seen that there is very little pebble diffusion, 
so fluctuations in the residence time are primarily due 
to hydrodynamic dispersion in the mean flow. We have 
also seen that the Kinematic Model gives a reasonable 
description of the mean flow profile in the conical fun- 
nel region, where most of the shear and hydrodynamic 
dispersion occur. Therefore, we can approximate the 
residence-time distribution by the distribution of times to 
travel along different streamlines of the mean fiow, start- 
ing from different radial positions, tq, at a given height 
zq. Below we will compare such predictions, based on 
our numerical solutions to the Kinematic Model, to our 
DEM simulations for the two reactor geometries. 



FIG. 10: (Color online) Horizontal profiles of porosity at dif- 
ferent heights in the 30° reactor geometry. 



B. An Analytical Formula 



in the container, as shown in Fig. |^b). As z increases 
from 30d to 60d, the interior of the packing goes from 
being disordered to having a strong radial ordering, cen- 
tered at around z — 12d. The reason for this ordering 
is due to the presence of the guide ring high in the con- 
tainer, which keeps the fuel and moderator pebbles sep- 
arate. The ring, placed at = 14. 5d in the container, 
creates radial crystallization, which can then propagate 
very far downward, since the packing is plug-like for most 



We can obtain a simple, exact formula for the 
residence-time distribution in a somewhat different ge- 
ometry using the Kinematic Model, as follows. The sim- 
ilarity solution to Eq. 10} for a wide, flat bottomed silo 
draining to a point orifice at z = is 



u(r, z) 
v{r, z) 



Qr 



-r^/46z 



Q 



-r' /Abz 



(5) 
(6) 



10 



where u and v are the radial (horizontal) and downward 
velocity components and Q is a constant proportional 
to the total flow rate through the orifice. (This is just 
the classical Green function for the diffusion equation in 
two dimensions, where z acts like "time".) A slightly 
more complicated solution is also possible for a parabolic 
silo, but let us focus on the simplest case of Eqs. Q- 
which is a good approximation for a wide parabolic 
funnel, where the velocity near the walls is small, i.e. 
R > y/ibzQ. A more detailed analysis is not appropriate 
here, since a simple analytical solution does not exist for 
the actual reactor geometry of a conical funnel attached 
to straight cylinder. 

For the flow field in Eqs. lO-®, the trajectory of a 
Lagrangian tracer particle along a streamline is given by 



dr 

'dt 

dz 

'dt 



ii(r, z), r{t = 0)=ro (7) 
-vir,z), z{t^O)^zo (8) 



Combining these equations and integrating, we find that 
the streamlines are parabolae, z/zq = (r/ro)^, and that 
the residence time for a pebble starting at (ro,Zo) is 



Toiro, zq) 



2Q 



(9) 



Now we consider pebbles that are uniformly dis- 
tributed at a height zq in a circular cross section of radius 
R in the flow field Eqs. (01-10 The probability distribu- 
tion for the residence time of those pebbles is 

p{t\zo,R) = [ S{t - To(ro,2:o)) ^^^°f'^° (10) 
Jo T^R^ 

for T < T„i„(zo) 

462:0/^^'^ for T^in <T < Tmax (H) 
iov T > Tjnax{zo,R) 



where 



T-o(0, Zq) 



h4 

2Q 



MR,zo)-—e 



(12) 
(13) 



Once again, this solution is strictly valid for an infinitely 
wide and tall silo draining to a point orifice, and it is 
roughly valid for a parabolic hmnel, z/ zq — (r/R)'^, as an 
approximation of a conical funnel in the actual reactor 
geometry. We can further approximate the effect of a 
nearly uniform flow of speed vq to describe the upper 
cylindrical region by simply adding (z — zq)/vq to the 
residence time for a starting point z > zq. 

Although this analysis is for a modified geometry, we 
will see that it captures the basic shape of the residence- 
time distributions from the DEM simulations in a simple 
formula Hll|) . The probability density is sharply peaked 
near the shortest residence time, Tmin, corresponding to 



pebbles near the central axis traveling the shortest dis- 
tance at the largest velocity. The longer distance and 
(more importantly) the smaller velocity at larger radial 
positions cause strong hydrodynamic dispersion, result- 
ing a fat-tailed residence-time density which decays like 

1/t, up to a cutoff Tmax- 



C. Simulation Results 

For the DEM reactor simulations, we calculate the dis- 
tribution of times it takes for particles to drop from sev- 
eral different values of zq, adding in a weighting factor to 
take into account that shorter residence times are pref- 
erentially observed in the data set. 

Since we are primarily interested in the radioactive 
burnup, we concentrate on the residence times for the 
fuel pebbles, but for comparison, we also report re- 
sults for the moderator pebbles. Figure [TTT a) shows the 
residence-time probability densities for pebbles starting 
at z = 40d, 55d, 70d to exit the container for the 30° 
reactor geometry. The distributions for the moderator 
pebbles are quite narrow, showing all particles exit over 
a short time window. In contrast, the distributions for 
the fuel pebbles exhibit fat tails, as expected qualita- 
tively from the Kinematic Model approximation Hll() for 
a parabolic geometry. A closer analysis of the data con- 
firms that the longest waiting times are associated with 
pebbles passing close to the walls, especially near the 
corner between the conical and cylindrical wall sections, 
although there are no completely stagnant regions. 

Figure mb) shows corresponding plots for the 60° re- 
actor geometry. In general, the residence-time densities 
have similar shapes as for the 30° geometry, but they 
are much narrower and exhibit a small secondary peak 
far into the tail. Examining movies shows that this extra 
peak is due to a boundary layer of particles, roughly one- 
pebble thick, touching the 60° conical wall sliding down 
at a speed lower than the nearby bulk. This extra source 
of hydrodynamic dispersion could not be easily captured 
by a continuum model for the mean flow. A simple way 
to eliminate it would be to replace add an outer annu- 
lus of moderator pebbles (controlled by another guide 
ring at the top) , which would flow more slowly along the 
walls, leaving the fuel pebbles in a more uniform flow 
with smaller fluctuations. Another possibility would be 
to reduce the wall friction, which makes the flow more 
uniform, as discussed in the following section. 

Figure E| investigates the accuracy of the Kinematic 
Model in predicting the DEM residence-time distribu- 
tion. The total residence-time distribution for both fuel 
and moderator pebbles to exit the reactor from z = 40d 
in the 30° geometry is shown, and is compared with two 
predictions from the Kinematic Model, one making use of 
the analytic formula (|ll|l . and one making use of the nu- 
merical solution of the velocity profile. We use of the 
value b = 2.5d and calibrate the total flow to match 
the total flow from the DEM data. Both the numeri- 
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FIG. 11: (Color online) Residence-time probability densities 
for the time it takes particles to drop from a specific height z 
out of the container, for the thirty degree reactor geometry (a) 
and sixty degree reactor geometry (b) for fuel pebbles (solid 
lines) and for moderator pebbles (dashed lines). 
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FIG. 12; (Color online) Comparison of the residence time 
distributions between DEM simulation, numerical solution of 
the Kinematic Model, and the analytic formula. 



cal solution and the analytic formula can roughly cap- 
ture the overall shape of the DEM distribution, although 
neither achieves a good quantitative agreement, particu- 
larly in the tails. Since the analytic formula assumes all 
streamlines are parabolic, it fails to take into account the 
slow-moving particles that stay close to the wall, and it 
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FIG. 13: (Color online) Distribution of times to pass through 
the entire container for fuel pebbles (solid lines) and moder- 
ator pebbles (dashed lines). 



therefore predicts a cut-off in the residence time distri- 
bution which is much shorter than some of the observed 
residence times in the DEM simulation. The numeri- 
cal solution of the Kinematic Model accounts for this 
and provides a better match, although it is clear that a 
model correctly accounting for the flow of pebbles near 
the container walls may be required in order to achieve 
high accuracy. 



D. Residence times for the entire container 

We also considered the distribution of times for the 
particles to pass through the entire container. While the 
flow in the upper part of the reactor is essentially plug- 
like, boundary effects near the container walls and on the 
guide ring can have an appreciable effect on the pebble 
residence times, which we study here. Since it takes a 
long time for particles to pass through the entire con- 
tainer we made use of the two extended data sets, con- 
sisting of 1,427 snapshots for the thirty degree geometry 
and 1,249 snapshots for the sixty degree geometry. 

Figure [TSI shows the time distributions for pebbles to 
pass through the entire container. Apart from a large 
positive time shift, the curves are similar in form to those 
in Fig. ^2 However, for both geometries, we see sec- 
ond small peaks in the distributions for the moderator 
pebbles, corresponding to a slow- moving boundary layer 
of pebbles touching the guide ring. The sixty degree 
curve for the fuel pebbles also exhibits several undula- 
tions corresponding to multiple layers of pebbles crys- 
tallized against the outer wall, each moving at different 
speeds. 



VII. WALL FRICTION 

The behavior of pebbles near the walls is of signiflcant 
interest to reactor design, and to look into this further. 
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FIG. 14; (Color online) Streamlines for the half-size, monodis- 
perse geometries with wall friction (left) and without wall fric- 
tion (right). Arrows are proportional to the velocity vectors 
in selected horizontal slices. 
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FIG. 15: (Color online) Comparison of velocity profiles for 
simulations with and without wall friction for two different 
heights. 
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FIG. 16: (Color online) Comparison of number density pro- 
files at a = 60d for simulations with and without wall friction. 



can see from figure 1161 that the number density profile 
is more peaked close to the wall. Figure ^1 also shows 
that the radial ordering created by the guide ring is also 
significantly enhanced. While this is due in part to the 
more plug-like flow allowing packing effects to propagate 
further down, it is also due to the frictionless guide ring 
initially creating radial ordering. Thus it may be possible 
to tune the material properties of the guide ring (or the 
roughness of its walls) to enhance or reduce the radial 
ordering effects. 

Removing wall friction also has the effect of increasing 
radial ordering effects near the wall. Perhaps most sur- 
prisingly, removing wall friction results in a significant 
alteration of the flow in the interior of the packing, as 
shown by the two velocity profiles in figure lT^ for z = 18c?. 
While both velocity profiles must converge upon the ori- 
fice, we see that the velocity profile for the /i^ = 0.7 
case is significantly more curved than that for /i^ ~ 0. 
This also has the effect of preferentially speeding up the 
relative flux of fuel pebbles: with wall friction, the fuel 
pebbles make up 71.5% of the total mass flux, but with- 
out wall friction, this increases to 74.7%. 



VIII. BIDISPERSITY 



we investigated the effect of wall friction by comparing 
two simulations runs in the half-size geometry, with wall 
friction coefficients ji^^ — and fiw = 0.7. All other 
aspects of the simulation, including the interparticle in- 
teractions, were kept the same. 

Figure 1151 shows a comparison of flow profiles for the 
two simulations at two different heights. We see that the 
/iji, — simulation results in a significantly larger flow 
speed, with a mass flow rate of 104mT^^, as opposed to 
59.6mT^^ for — 0.7. As would be expected, removing 
wall friction also removes the boundary layer of slower 
velocities at the wall, creating an almost perfectly uni- 
form velocity profile high in the reactor. This also has 
the effect of increasing radial ordering effects, and we 



A. The Bidisperse PBR Concept 

The two-pebble design of MPBR with a dynamic cen- 
tral moderator column has various advantages over a 
solid graphite central column (as in the revised PBMR 
design). For example, it flattens the neutron flux profile, 
while preserving a very simple core vessel without any 
internal structures, which would be subjected to extreme 
radiation and would complicate the granular flow. It also 
allows the widths of the moderator column and fuel an- 
nulus to be set "on the fly" during reactor operation, 
simply by adjusting the guide ring at the top. 

A drawback of the dynamic moderator column, how- 
ever, is its porosity, which allows the passage of the 
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FIG. 17: Schematic diagram of the pebble flow in a bidisperse 
MPBR design. 



helium-gas coolant, at the highest velocity (along the cen- 
tral axis). To improve the thermal efficiency and power 
output, it would be preferable to focus the gas flow on 
the fuel annulus and the interface with the moderator 
column, where the most heat is generated. This is auto- 
matically achieved with a solid graphite central column, 
but there is a very simple way to shape the gas flow in a 
similar way with a dynamic column, while preserving its 
unique advantages. 

The idea is to make the graphite moderator pebbles in 
the central column smaller than the fuel pebbles in the 
outer annulus, as shown in Fig. 1171 (This also helps with 
sorting of fuel and moderator pebbles as they exit the 
core.) In standard continuum models of flow in porous 
media [33 . l35l l3^ , the permeability of the packing scales 
with the square of the pebble diameter (or pore size), 
so reducing the diameter of the moderator pebbles can 
greatly reduce the gas flow (e.g. by a factor of four for 
half-diameter pebbles). This argument holds everywhere 
that the packing is statistically the same, in the monodis- 
perse packings of the fuel annulus and the moderator 
column, which have the same porosity. At the interface 
between the two regions, we have seen in Figures Eland 
^]that the porosity is enhanced for a monodisperse core 
due to the guide ring, although a bidisperse interface will 
have different structure. 

In any case, it is clear that the bidisperse core will fo- 
cus the coolant flow away from the moderator column 
and onto the fuel annulus, as shown in Figure ^| In 
most PBR designs, high-pressure helium gas is intro- 
duced from a reservoir above the core, through holes in 
the graphite bricks which make up the core vessel. The 



graphite bricks enclosing the core 




FIG. 18: Sketches of the helium-gas coolant system (top) 
and the flows a basic PBR core (left), a monodisperse MPBR 
where gas is introduced only in the fuel annulus outside the 
guide ring (center), and a bidisperse MPBR where the dy- 
namic moderator column has much lower permeability due to 
smaller pebble size (right). 



gas then flows through the core and exits through holes in 
the graphite bricks in the conical funnel to another reser- 
voir at "very high" temperature (w 950° C). In MPBR, 
the gas can be introduced only outside the guide ring, 
which focuses the gas flow on the fuel annulus down to a 
distance comparable to the radius of the guide ring. With 
a significant reduction in permeability of the central col- 
umn in the bidisperse core, the gas flow can be focused 
almost entirely on the fuel annulus and the interfacial 
region (where the heat generation is maximal). 



B. Simulation Results 

The only question regarding the feasibility of the bidis- 
perse core is the stability of the central column over time 
and the possibility of enhanced diffusion of the small 
moderator pebbles into the annulus of larger fuel pebbles. 
In other systems, such as rotating drums Isol Isif . 



14 




FIG. 19: (Color online) Snapshots of vertical cross-sections 
for the bidisperse simulations. From left to right, the moder- 
ator pebbles have diameters Id, 0.8d, and 0.5d while the fuel 
pebbles are of constant size Id. 

vibrated buckets [s^ Is^ , and draining silos bidis- 
perse granular materials display a tendency to segregate 
(rather than mix) during dynamics, but there is currently 
no general theory which could be applied to our reactor 
geometry. Therefore, our DEM simulations provide a 
useful means to address this important question. 

Figure [TOI shows snapshots of vertical cross sections for 
the three different bidisperse simulations that were run 
in the half-size geometry. As shown in the diagram, the 
central column remains stable and coherent in all three 
cases, and very little mixing between the two types of 
pebbles is visible. Figure 1^ shows a comparison of the 
velocity profiles from the three simulations for two dif- 
ferent heights. It is reassuring to see that the bidisperse 
simulations do not significantly differ from the monodis- 
perse simulation, although we do see a slightly higher 
overall flow rate in the bidisperse systems: we see total 
mass flow rates of 59. Gmr"^, 60.8mT~^, and 65.0mr~^ 
for the monodisperse, 0.8:1, and 0.5:1 simulations respec- 
tively. 

The velocity profiles are slightly more curved in the 
bidisperse central core; this is particularly apparent in 
the 0.5:1 simulation. This leads to a small cusp in the 
velocity profile near the interface between the two types 




2 4 6 8 10 12 14 
Radial distance (d) 

FIG. 20: (Color online) Comparison of velocity profiles for 
the three bidisperse simulations. The three flatter curves are 
calculated at z = 30d in the plug-like flow region while the 
other three were taken at z — 22d in the parabolic flow region. 
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FIG. 21: (Color online) Comparison of particle diffusion for 
the three bidisperse simulations. 



of particles which may lead to adverse mixing effects. The 
faster flow also leads to a significantly larger turnaround 
of the moderator pebbles. In the monodisperse system, 
the moderator pebbles comprise 28.5% of the total mass 
flux, but this is increased to 31.7% in the 0.8:1 bidisperse 
simulation, and 42.6% in the 0.5:1 bidisperse simulation. 

To investigate the amount of mixing of the central col- 
umn, we used a technique similar to that described in 
section Hvl At z = llOo? all moderator particles with 
r > 8d are marked, and their radial diffusion is then cal- 
culated as a function of z. The results are shown in figure 
1211 in the cylindrical section of the packing, there is very 
little difference between the three simulations, but in the 
area of convergent flow, we see that bidispersity leads to 
significantly more mixing. However, even for the 0.5:1 
simulation, the scale of diffusion is still smaller than a 
single particle diameter, and essentially the central col- 
umn remains stable. 

Due to computational limitations, we were unable to 
investigate smaller size ratios in the reactor geometries, 
so we carried out simulations in a smaller container with 
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a 0.3:1 size ratio and found dramatically different be- 
havior: During drainage, the central column became un- 
stable, and the small particles penetrated many particle 
diameters into the packing of larger particles. We expect 
that there is a fundamental crossover in behavior sim- 
ply due to geometry of amorphous packings, when the 
moderator pebbles become small enough to pass through 
the gaps between the densely packed fuel pebbles. An 
in-depth study of this phenomenon remains a subject of 
future work. For now, we can safely recommend a di- 
ameter ratio of 0.5:1, which reduces the dynamic central 
column's permeability by a factor of four without intro- 
ducing any significant diffusion of moderator pebbles into 
the fuel annulus. 



IX. CONCLUSIONS 
A. Pebble-Bed Reactor Core Design 

Using DEM simulations, we have analyzed many as- 
pects of granular flow in pebble-bed reactor cores of di- 
rect relevance for design and testing. We close by sum- 
marizing some key conclusions. 

The mean flow profile exhibits a smooth transition 
from a nearly uniform plug flow in the upper cylindri- 
cal region to a nonuniform, converging flow in the lower 
funnel region, consistent with recent experiments 9, 37]. 
There are no stagnant regions in the 30° and 60° con- 
ical funnels considered in this study, although the flow 
is slower near the corner at the top of the funnel, espe- 
cially in the former case. Moreover, the wider 30° funnel 
has a boundary mono-layer of slower pebbles partially 
crystallized on the wall. 

The only available continuum theory for such flows, 
the simple Kinematic Model [13, llS, U2, ; gives a rea- 
sonable qualitative picture of the flow profiles, although 
it cannot capture discrete boundary-layer effects. As in 
other experiments on similar geometries P], the Kine- 
matic Model does not quantitatively predict the depen- 
dence of the flow profile on geometry. We suggest that 
it be used to get a rough sense of the flow profile for a 
given core geometry prior to (much more computation- 
ally expensive) DEM simulations and/or experiments. 

We have quantified the degree of pebble mixing in the 
core. Although there is some horizontal diffusion in the 
funnel region, pebbles depart from the streamlines of the 
mean flow by less than one pebble diameter prior to ex- 
iting the core. 

We have demonstrated that the "mixing layer" be- 
tween the central moderator column and the outer fuel 
annulus, which appears in prior models "3^, can be re- 
duced to the thickness of one pebble diameter by sepa- 
rating moderator and fuel pebbles with a guide ring at 
the ceiling (to eliminate mixing by surface avalanches), 
consistent with experiments on MPBR models [s^- We 
conclude that the dynamic central column of moderator 
pebbles is a sound concept, which should not concern 



regulators. 

We have constructed Voronoi tessellations of our flow- 
ing packings to measure the profile of volume fraction 
(or porosity) and found some unexpected features which 
would affect coolant gas flow through the core. The bulk 
of the core, in the plug-flow region of the upper cylinder, 
has a volume fraction near the jamming point (63%), 
but there is a sharp transition to less dense packings 
(55 — 60%) in the funnel region, due to shear dilation. We 
also observe lower volume fractions in this range at the 
moderator/fuel interface in the upper cylinder, below the 
guide ring, and lower volume fractions (50 — 55%) against 
the walls. These narrow regions of increased porosity 
(and thus, increased permeability) would allow faster he- 
lium gas flow. 

We have also studied local ordering in the flowing pack- 
ings and flnd evidence for partial crystallization within 
several pebble diameters of the walls, consistent with pre- 
vious experiments 32. 33J and simulations |40l . Such 
ordering on the walls of the guide ring, then advected 
down through the core, is responsible for the increased 
porosity of the moderator/fuel interface. 

We have varied the wall friction in our DEM simula- 
tions and observe that it can affect the mean flow, even 
deep into the bulk. Reducing the wall friction increases 
radial ordering near the walls and makes the flow proflle 
more uniform. 

Since diffusion is minimal, the probability distribution 
of pebble residence times is dominated by advection in 
the mean flow. Therefore, we have made predictions us- 
ing the Kinematic Model, numerically for the conical- 
funnel reactor geometries, and analytically for a wide 
parabolic funnel. The model predicts a fat-tailed (~ 1 /t) 
decay of the residence-time density due to hydrodynamic 
dispersion in the funnel region. 

Our DEM simulations predict that the 60° conical fun- 
nel results in a narrower residence-time distribution than 
the 30° funnel, which has more hydrodynamic dispersion. 
The steeper 60° funnel also exhibits a boundary layer of 
slower, partially crystallized pebbles near the wall which 
lead to an anomalous bump far in the tail of residence- 
time distribution. These results have important implica- 
tions for non-uniformity in the burnup of fuel pebbles. 

We have introduced the concept of a bi-disperse core 
with smaller moderator pebbles in the dynamic central 
column than in the outer fuel annulus, in order to fo- 
cus the helium gas flow on the fuel. Our DEM simula- 
tions demonstrate that there is negligible pebble mixing 
at the interface for diameter ratios as small as 0.5:1, for 
which the permeability of the moderator column is re- 
duced by a factor of four. We conclude that the bidis- 
perse MPBR design is sound and will produce a stable 
moderator-pebble column of greatly reduced gas perme- 
ability. 

A natural next step would be to combine our full-scale 
DEM model for the pebble flow with existin g co mputa- 
tional approaches to reactor core physics 'SS'/S^, which 
rely on pebble flow as an empirical input. More accurate 
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studies of gas flow in the core could also be done, start- 
ing from our complete pebble packings, or the average 
quantities such as the porosity. With such computational 
tools, one should be able to reliably test and develop new 
reactor designs. 



B. Basic Physics of Dense Granular Flow 

We have noted a number of favorable comparisons be- 
tween our simulations and experiments in similar geome- 
tries 0, li^, 113, 113 , which provides further validation of 
the Discrete-Element Method as a realistic means of sim- 
ulating granular materials. As such, it is interesting to 
consider various implications of our results for the the- 
ories of dense granular flow, since the simulations probe 
the system at a level of detail not easily attained in ex- 
periments. 

Our conclusions about the Kinematic Model are simi- 
lar to those of a recent experimental study 9] : The model 
describes the basic shape of the flow field in the converg- 
ing region, but fails to predict the nearly uniform plug 
flow in the upper region with vertical walls or the precise 
dependence on the funnel geometry. It also cannot de- 
scribe boundary-layers due to partial crystallization near 
walls or incorporate wall friction, which we have shown 
to influence the entire flow profile. 

On the other hand, there is no other continuum 
model available for dense silo drainage, except for Mohr- 
Coulomb plasticity solutions for special 2d geometries, 
such as a straight 2d wedge without any corners [T^ . 
so it is worth trying to understand the relative success 
of the Kinematic Model for our 3d reactor geometries 
and how it might be improved. A cooperative micro- 
scopic mechanism for random-packing dynamics, based 
on "spots" of diffusing free volume, has recently been 
proposed, which yields the mean flow of the Kinematic 
Model as the special case of independent spot random 
walks with uniform upward drift from the orifice (due 
to gravity) [l7|- Under the same assumptions, the Spot 
Model has also been shown to produce rather realistic 
simulations of flowing packings in wide silos (compared 
to DEM simulations) [l8|. where the Kinematic Model 
is known to perform well [T^ [Tj, [l^ [ig ■ This suggests 
that some modification of the spot dynamics, such as 
spot interactions and/or nonuniform properties coupled 
to mechanical stresses, and an associated modification 
of the Kinematic Model in the continuum limit, may be 
possible to better describe general situations. 

From a fundamental point of view, perhaps the most 
interesting result is the profile of Voronoi volume fraction 
(or porosity) in our fiowing random packings in Figure 
Although the mean velocity in Figure |21 shows a fairly 
smooth transition from the upper plug flow to the lower 
converging flow, the volume fraction reveals a sharp tran- 
sition (at the scale of 1 — 3 particles) from nearly jammed 
"solid" material in the upper region (63%) to dilated, 
sheared "liquid" material in the lower region (57-60%). 



The transition line emanated from the corners between 
the upper cylinder and the conical funnel. We are not 
aware of any theory to predict the shape (or existence) 
of this line, although it is reminiscent of a "shock" in the 
hyperbolic equations of 2d Mohr-Coulomb plasticity [T^ . 

Our measurements of diffusion and mixing provide 
some insights into statistical fluctuations far from equilib- 
rium. Consistent with the experiments in wide quasi-2d 
silos 10 , we find that diffusion is well described geomet- 
rically as a function of the distance dropped, not time 
(as in the case of thermal molecular diffusion). As a 
clear demonstration, there is essentially no diffusion as 
pebbles pass through the upper core, until they cross the 
transition to the funnel region, where the diffusion re- 
mains small (at the scale of one pebble diameter) and 
cooperative in nature. The behavior in the funnel is con- 
sistent with the basic Spot Model 0, but a substantial 
generalization would be needed to describe the transition 
to the upper region of solid-like plu g fiow, perhaps using 
concepts from plasticity theory [26j. 

We view silo drainage as a fundamental unsolved prob- 
lem, at least as interesting and relevant for applications 
as Couette shear cells, which have been received much 
more attention in physics. The challenge will be to find 
a single theory which can describe both shear cells and 
silo drainage. Our results for pebble-bed reactor geome- 
tries may provide some useful clues. 
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APPENDIX A: NUMERICAL SOLUTION OF 
THE KINEMATIC MODEL 

In the Kinematic Model for drainage the vertical down- 
ward velocity u in the container is assumed to follow a 
diffusion equation of the form 



where V^^ is the horizontal Laplacian. By exploiting the 
axial symmetry, v can be treated as a function of z and 
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r only. In cylindrical coordinates the Laplacian is 



dv 
'd'z 



,1 d { dv 
0-^ r 



r dr 



dr 
1 dv 
r dr 



To numerically solve this partial differential equation, 
we make use of the implicit Crank-Nicholson integration 
scheme. We write = ^(j AA, nA?/), and solve in the 
range j = 0, 1, . . . , iV where N = AA~^. Away from the 
end points, the Crank-Nicholson scheme tells us that 



The radial velocity component is given by 

, dv 
dr 

and by enforcing that the velocity field at the wall must 
be tangential to the wall, we can obtain boundary con- 
ditions for solving v. 

To solve the above equation in a cylinder is straight- 
forward, since we can make use of a rectangular grid. 
The boundary condition reduces to Vr = at the wall. 
However, to solve this equation in the reactor geometry, 
we must also consider the complication of the radius of 
the wall, R, being a function of z. To ensure accurate 
resolution in the numerical solution of v at the wall, we 
introduce a new coordinate A = r/R{z),r] = z, which 
then allows us to solve for u over the range < A < 1. 
Under this change of variables, the partial derivatives 
transform according to 

d_ _ 1 d 
dr Riv) 9^ 

d_ _ d XR'jr]) d 
dz dy R{v) 9X 

In the transformed coordinates 

b 



R^Vr, 



A 



v\ + bvxx + XRR'vx- 



To ensure differentiability at r = 0, we use the boundary 
condition 



A77 



2AA2i?2 ^ J+i 



,,"+1 



where all references to R and R' are evaluated at 77 = 
A?7(j -I- i). If j = 0, then by reference to equation lAll 
we find that 



,,"+1 _ „" 



."+1 _ ,,"+1 _1_ _ n,^\ 
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Similarly, for j = A^, by reference to eauation lA2l we see 
that effectively 



v'f^R'R 



2AA 



and hence 



7;"+^ - V 



N 



A?] 



v';,+' + v 



AA2i?2 

'{2N+1)R' i?'2 
2R ^ "26" 
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dv 
dX 



= 0, 



(Al) 



A=0 



and by ensuring zero normal velocity at the wall we find 
that 



dv 
dX 



vR'R 



(A2) 



If we write v" = (uq . . . then the above nu- 

merical scheme can be written in the form S'v"^^ = Tv" 
where S and T are tridiagonal matrices; this system can 
be efficiently solved by recursion in 0{N) time. The 
above scheme was implemented in CH — h, and gives ex- 
tremely satisfactory results, even with a relatively small 
number of gridpoints. 
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